【题解】 [国家集训队]Crash的数字表格莫比乌斯反演 luoguP1829

吐槽一下Typora这个编辑器：码了一上午的题解，居然突然卡机，并且自动关掉了，然后重新打开，发现保存的也没了。然后弹出一个“Typora意外关闭”的窗口，真想一拳上去。只好重新自己码了……。(以上是吐槽，请不要在意)

算了算了，重新写吧。所以你看到的这是第二份稿子。

仍然上莫比乌斯反演。

众所周知：

$lcm(i,j)=\frac{ij}{gcd(i,j)}$

那么我们将这个带进原式：

$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) = \sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}$

我们枚举 $gcd(i,j)$ 的值：

$\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]\frac{ij}{d}$ $=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]ijd$ $=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]ij$

设：

$f(x)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=x]ij$ $g(x)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[x|gcd(i,j)]ij$

可得：

$g(x)=\sum_{x|d}f(d)$

考虑怎么计算 $g(x)$ ：

$g(x)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[x|gcd(i,j)]ij$ $g(x)=\sum_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dx}\rfloor}ij \cdot x^2$ $g(x)=x^2\sum_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dx}\rfloor}ij$

这个显然是可以 $O(1)$ 算出的。

继续：

$ans=\sum_{d=1}^{n}d\cdot f(1)$ $f(1)=\sum_{d=1}^{\lfloor\frac{n}{}\rfloor}\mu(d)g(d)$

这个时候的复杂度只是 $O(n^2)$ ，继续优化。

将 $ans$ 写出：

$ans=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)=1]ij$

发现后面可以整除分块！

继续，将 $f(1)$ 写出：

$f(1)=\sum_{i=1}^{n}\mu(i)g(i)$ $f(1)=\sum_{i=1}^{n}\mu(i)i^2\sum_{a=1}^{\lfloor\frac{n}{di}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{di}\rfloor}ab$

将 $g$ 拆开后，我们可以发现后面又可以整除分块！

那么现在就是 $O(n)$ 了，可以过。

Code-$O(n)$

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>

#define ll long long
#define MOD 20101009
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))

const int N=1e7+2;
const int inf=1e9+9;

int mui[N],sum[N];
int vis[N],prime[N],cnt;

inline void _pre_mui(int n){
    mui[1]=1;
    for(int i=2;i<=n;++i){
        if(!vis[i])prime[++cnt]=i,mui[i]=-1;
        for(int j=1;j<=cnt;++j){
            if(i*prime[j]>n)break;
            vis[i*prime[j]]=1;
            if(!(i%prime[j])){mui[i*prime[j]]=0;break;}
            mui[id*prime[j]]=-mui[i];
        }
    }return;
}

inline int solve(int n,int m){
    ll ans=0;
    for(int l=1,r=0;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        int res=1ll*(1ll*(n/l)*(n/l+1)/2%MOD)*(1ll*(m/l)*(m/l+1)/2%MOD)%MOD;
        ans+=1ll*(sum[r]-sum[l-1])%MOD*res%MOD;
        ans%=MOD;
    }return (ans+MOD)%MOD;
}

int main(){
    int n,m,ans=0;
    scanf("%d%d",&n,&m);
    if(n>m)n^=m^=n^=m;
    _pre_mui(n);
    for(int i=1;i<=n;++i)
        sum[i]=(sum[i-1]+1ll*i*i%MOD*mui[i]%MOD+MOD)%MOD;
    for(int l=1,r=0;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        int res=1ll*(l+r)*(r-l+1)/2%MOD;
        ans=(ans+1ll*solve(n/l,m/l)*res%MOD)%MOD;
    }
    printf("%d\n",ans);
    return 0;
}

毒瘤出题人不会放过我们，这个毒瘤更改了数据： $n,m \leq 10^{10}$ 。

~~“哇咔咔咔卡掉你们的 O(n) !”~~

~~真想一拳上去这个智障。~~

$O(\sqrt{n})$ 是过的了的，故考虑向着这个方向前进。

继续推式子：

$ans=\sum_{d=1}^{n}d\sum_{i=1}^{n}\mu(i)i^2\sum_{a=1}^{\lfloor\frac{n}{di}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{di}\rfloor}ab$

看见这个 $di$ 了吗？我们令 $T=di$ ，然后将 $T$ 扔到前面去枚举一下。

然后就是后面的两个 $\sum$ ，这玩意跟 $i,d$ 没关系，一起扔到前面去。

于是就变成了：

$ans=\sum_{T=1}^{n}\sum_{a=1}^{\lfloor\frac{n}{T}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{T}\rfloor}ab\sum_{i|T}i^2\frac{T}{i}\mu(i)$

？？？？？？？

首先前面的这段是没有问题的对吧，那么后面的呢？

后面的原来是不是：

$\sum_{d=1}^{n}d\sum_{i=1}^{n}\mu(i)i^2$

那么……现在我们枚举的 $T$ 是 $i \cdot d$ ，我们枚举了一下可能的 $i$ ，成为 $i$ 的必备条件肯定是能被 $T$ 整除对吧？那么这个时候的 $d$ 呢？很显然是 $\frac{T}{i}$ 对吧？

所以啊……就是这么写了。

但是这样写有什么用啊 $QwQ$

首先，我们来看，$i^2$ 是个什么鬼？我们设一个函数 $Q(i)=i^2$ ，于是我们可以发现，这个东西是个完全积性函数，然后看 $\frac{T}{i}$ ，显然是 $id(\frac{T}{i})$ ，也是完全积性函数。于是两个完全积性函数用狄利克雷卷积卷起来，它们的狄利克雷卷积是一定可以筛出来的。后面的 $\mu$ 是积性函数，然后呢？后面的那一坨都可以筛出来！

于是美滋滋。

我们设 $sum[k]$ 表示当 $T$ 为 $k$ 的时候后面那一坨的值。

那么现在分三种情况：

$k$ 是质数，这下子后面的 $i$ 只能是 $1$ 和 $k$ ，$1$ 的时候的值就是 $k$ ，$k$ 的时候的值是 $k^2\cdot \frac{k}{k}\cdot \mu(k)$ ，很显然这个时候的 $\mu(k)$ 的值是 $-1$ ，于是这个时候的值是 $-k^2$ ，那么这个时候 $sum[k]$ 的值是 $k-k^2$ 。
$\mu(k)$ 为 $0$ ，这下子的话就肯定有一个 $j$ ，使得 $k$ 可以整除 $j^2$ ，这个时候假设就只能整除 $j^2$ ，也就是说 $\mu(k/j)$ 的值非 $0$ 。那么我们看看，在 $sum$ 所计算的式子中，只有 $T$ 的因子对 $T$ 产生贡献。考虑 $k/j$ 到 $k$ 多了什么因子。这个时候多的因子有两类，一类是包含了 $j^2$ 的，一类是只包含了 $j$ 的。第二类的可以先不管，因为之前 $k/j$ 中有了一个 $j$ ，这类因子的贡献已经算过了。那么对于第一类因子，因为包含了 $j^2$ ，所以 $\mu$ 值为 $0$ ，对答案没有任何贡献。

那么这个时候对答案有贡献的还是 $k/j$ 的因子，乘上一个 $j$ 后没有更多的对答案造成贡献的因子。

但是我们发现上限 $T$ 变了，增大了 $j$ 倍，对于原来的每份贡献的值也增大了 $j$ 倍。由于没有其他的贡献，$k$ 的所有的贡献都来自 $k/j$ ，那么直接转移就好。

所以是$sum[k]=sum[k/j]\cdot j$
对于剩下的情况，我们发现，这个可以直接转移了。当我们枚举 $k$ 的时候，考虑怎么用 $k$ 来转移 $k\cdot j$ ，这个时候 $j$ 是质数，并且 $k$ 中不包含 $j$ ，也就是说 $k$ 与 $j$ 互质。于是根据积性函数的性质，$sum[k]=sum[k/j]\cdot sum[j]$ 就好。

于是这个时候前面再整出分块一下，复杂度 $O(\sqrt{n})$ 。

~~听说有人被卡住 $O(n)$ 后没有推式子了，直接上了个杜教筛，这人一看就是杜教士了，并且也说明不珂学的上杜教筛是布星的~~

Code-$O(\sqrt{n})$

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>
 
#define ll long long
#define MOD 20101009
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
 
const int N=1e7+3;
const int inf=1e9+9;
 
int sum[N],vis[N],prime[N],cnt;
 
inline void _pre_mui_sum(){
    vis[1]=sum[1]=1;
    for(int i=2;i<=N;++i){
        if(!vis[i])prime[++cnt]=i,sum[i]=(i-1ll*i*i%MOD+MOD)%MOD;
        for(int j=1;j<=cnt;++j){
            if(i*prime[j]>=N)break;
            vis[i*prime[j]]=1;
            if(!(i%prime[j])){sum[i*prime[j]]=1ll*sum[i]*prime[j]%MOD;break;}
            else sum[i*prime[j]]=1ll*sum[i]*sum[prime[j]]%MOD;
        }   
    }
    for(int i=1;i<=N;++i)
        sum[i]=(sum[i-1]+sum[i])%MOD;
}
 
int main(){
    int n,m;
    _pre_mui_sum();
    scanf("%d%d",&n,&m);
    if(n>m)std::swap(n,m);
    long long ans=0;
    for(int l=1,r=0;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        int res=(1ll*(1+n/l)*(n/l)/2%MOD)*(1ll*(1+m/l)*(m/l)/2%MOD)%MOD;
        ans+=1ll*(sum[r]-sum[l-1]+MOD)%MOD*res%MOD;
        ans%=MOD;
    }
    printf("%lld\n",(ans+MOD)%MOD);
    return 0;
}

本文标题:【题解】 [国家集训队]Crash的数字表格莫比乌斯反演 luoguP1829

文章作者:Qiuly

发布时间:2019年03月01日 - 00:00

最后更新:2019年03月29日 - 13:52

原始链接:http://qiulyblog.github.io/2019/03/01/[题解]luoguP1829/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际转载请保留原文链接及作者。